QR_Householder_decomposition Subroutine

public subroutine QR_Householder_decomposition(A, Q, R)

QR decomposition using Householder method

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), DIMENSION(:, :) :: A
real(kind=dp), intent(out), DIMENSION(SIZE(A, 1) ,SIZE(A, 2)) :: Q
real(kind=dp), intent(out), DIMENSION(SIZE(A, 1) ,SIZE(A, 2)) :: R

Calls

proc~~qr_householder_decomposition~~CallsGraph proc~qr_householder_decomposition QR_Householder_decomposition proc~identity_n Identity_n proc~qr_householder_decomposition->proc~identity_n

Called by

proc~~qr_householder_decomposition~~CalledByGraph proc~qr_householder_decomposition QR_Householder_decomposition proc~qr_decomposition QR_decomposition proc~qr_decomposition->proc~qr_householder_decomposition proc~eigen Eigen proc~eigen->proc~qr_decomposition proc~is_spd is_SPD proc~is_spd->proc~eigen

Source Code

    SUBROUTINE QR_Householder_decomposition(A, Q, R)
        REAL(dp), DIMENSION(:, :), INTENT(IN) :: A
        REAL(dp), DIMENSION(SIZE(A, 1) ,SIZE(A, 2)), INTENT(OUT) :: Q, R
        REAL(dp), DIMENSION(SIZE(A, 1) ,SIZE(A, 2))::Id, H, v_mat_tmp
        REAL(dp), DIMENSION(SIZE(A, 1)) :: v, u, x
        INTEGER :: N, i, j, k
        REAL(dp) :: alpha, w, signe, norm_u

        N= SIZE(A, 1)

        R = A

        Id = Identity_n(N)
        Q = Identity_n(N)

        DO k = 1, N

            x = 0.d0
            u = 0.d0
            v = 0.d0
            v_mat_tmp = 0.d0
            x(k:N) = R(K:N, K)

            alpha = NORM2(R(k:N, k))

            signe = - SIGN(alpha,x(k))
            u(k:N) =  x(k:N) - signe * Id(k:N, k)

            norm_u = NORM2(u)
            IF (norm_u < epsi) CYCLE
            v(k:N) = u(k:N) / norm_u

            w = 1.d0
            DO i = k, N
                DO j = k, N
                    v_mat_tmp(i, j) = v(i) * v(j)
                END DO
            END DO

            H = Id
            H(k:N, k:N) = Id(k:N, k:N) - (1.d0 + w) * v_mat_tmp(k:N, k:N)

            Q = MATMUL(Q, H)
        
            R(k:N, k:N) = MATMUL(H(k:N, k:N), R(k:N, k:N))
    
        END DO
    END SUBROUTINE QR_Householder_decomposition